We are presenting an Internal Linear Combination (ILC) CMB map, in which the foreground is reduced through harmonic variance minimization. We have derived our method by converting a general form of pixel-space approach into spherical harmonic space, maintaining full correspondence. By working in spherical harmonic space, spatial variability of linear weights is incorporated in a self-contained manner and our linear weights are continuous functions of position over the entire sky. The full correspondence to pixel-space approach enables straightforward physical interpretation on our approach. In variance minimization of a linear combination map, the existence of a cross term between residual foregrounds and CMB makes the linear combination of minimum variance differ from that of minimum foreground. We have developed an iterative foreground reduction method, where perturbative correction is made for the cross term. Our CMB map derived from the WMAP data is in better agreement with the WMAP best-fit $\Lambda$CDM model than the WMAP team's Internal Linear Combination map. We find that our method's capacity to clean foreground is limited by the availability of enough spherical harmonic coefficients of good Signal-to-Noise Ratio (SNR).